A model for gelation with explicit solvent effects: Structure and dynamics 
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We study a two-component model for gelation consisting of /-functional monomers (the gel) and 
inert particles (the solvent). After equilibration as a simple liquid, the gel particles are gradually 
crosslinked to each other until the desired number of crosslinks has been attained. At a critical 
crosslink density the largest gel cluster percolates and an amorphous solid forms. This percolation 
process is different from ordinary lattice or continuum percolation of a single species in the sense 
that the critical exponents are new. As the crosslink density p approaches its critical value pc, the 
shear viscosity diverges: rj(p) ~ (pc—p)~" with s a nonuniversal concentration-dependent exponent. 

PACS numbers: 61.43.Hv,66.20.-|-d,83.10.Rs,83.60.Bc 
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I. INTRODUCTION 

It is generally accepted that percolation is an essen- 
tial aspect of gelation or vulcanization — it is doubtful 
that even in a highly entangled melt of long polymers a 
nonzero value of the static shear modulus could exist in 
the absence of an infinite connected network. However, 
percolation has usually been studied in rather special lim- 
its. Site and bond percolation of a single species on reg- 
ular lattices are very well characterized and off-lattice 
percolation seems to present no new features , at least 
insofar as critical behavior is concerned. More closely re- 
lated to real gels are the so-called correlated percolation 
models where the distribution of crosslinks is drawn from 
a Boltzmann distribution appropriate for a nearest neigh- 
bor lattice gas Except at special points in the phase 
diagram these models are also in the universality class of 
the simple percolation problem. In our previous work on 
transport properties near the gel point ||^, we have also 
used a simple one-species percolation process to produce 
the incipient gel. We found that the shear viscosity di- 
verges as the percolation concentration pc is approached 
according to rj{p) ~ {pc — p)~^ with s ~ 0.7. This value 
of the exponent s is in excellent agreement with a pre- 
diction of de Gennes based on a superconductor-normal 
conductor analogy Q and with recent analytical work on 
a Rouse model . It is also reasonably close to some ex- 
perimental results for s |^ but quite different from that 
produced by another set of experiments 1.1 < s < 1.3 
§. Thus it seems reasonable to ask if different versions 
of the crosslinking process might produce significantly 
different cluster size distributions from percolation and, 
consequently, different rheological properties. 

Gelation often occurs in the presence of a solvent and 
over some period of time rather than instantaneously, as 
in the usual percolation models. To simulate this fea- 
ture, we have considered a two-species model consisting 
of a fraction c of /-functional particles that are eligible 
to bond irreversibly to others of the same kind. The re- 
maining particles are inert and function as a background 



liquid through which the gel particles and clusters diffuse. 
Crosslinking occurs in stages: the equations of motion of 
all the particles are integrated forward for a fixed number 
of time steps between crosslinking attempts and this pro- 
cess is continued until the desired number of crosslinks is 
attained. At a critical concentration of crosslinks, Pc, (in 
the thermodynamic limit) the largest cluster percolates 
and an amorphous solid forms. For this process one can 
calculate the usual static or geometrical quantities used 
to characterize percolating systems, e.g., the fraction of 
particles on the 'infinite cluster', Poo(p) ^ {p — Pc)^ , the 
mean mass of finite clusters, S{p) ~ \p — Pc\^^ , the frac- 
tion of samples percolating f{p), the cluster size distribu- 
tion n{m,p) = m~'^ (j){m\p — Pc\^^'^) where m is the mass 
of a cluster and the radius of gyration Rg{m) ~ m}/^ 
where D is the fractal dimension of the clusters. For 
simple percolation processes, r « 2.18, tr ~ 0.45 and 
these two exponents determine the others through scal- 
ing relations 0|. Here we find, at least for small c, that 
the cluster size distribution, even at pc, is not well de- 
scribed by a simple power law. However, the other static 
quantities listed above do display power law behavior 
near pc and a standard finite-size scaling analysis pro- 
vides a very good collapse of our data. Moreover, the 
hyperscaling relation 2/3 -I- 7 = dv, where d = i is the 
dimensionality and v the correlation length exponent, is 
satisfied. This suggests that this percolation transition is 
fundamentally describable in terms of a fixed point with 
two (at least) relevant scaling fields. As the percolation 
point is approached from below, the shear viscosity di- 
verges according to T]{p) ^ {pc —p)~'^. In contrast to our 
previous work on a model without solvent, we find val- 
ues of s in the range 0.3 < s(c) < 0.45 as compared with 
s ~ 0.7. These results suggest that the critical behavior 
of transport coefficients of systems close to the gel point 
is nonuniversal. 

The structure of this article is as follows. In section II 
we describe the present model and simulation procedures 
in more detail. The geometric properties of the system 
are discussed in section HI and the data on the shear 
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viscosity are presented in section [V. We conclude with 
a brief summary and discussion in section 0. 



II. THE MODEL 

We consider a system of N particles in three dimen- 
sions, all of which interact with each other through the 
soft-sphere potential V{rij) = e((To/rjj)^^ for Vij < I.Sctq 
1^ with (To = 1 and ksT/e = 1. We simulated sys- 
tems at a volume fraction $ — ttctqA^/GV^ = 0.4 which 
is well below the liquid-solid coexistence density. In the 
absence of any other interactions, this system would be 
a simple three-dimensional liquid. We initially place the 
particles on a simple cubic lattice that fills the computa- 
tional box. We then randomly select Ngei — cN particles 
to be the gel forming component. After equilibration 
of the system with Brownian dynamics, with periodic 
boundary conditions, for 10000 time steps we begin the 
crosslinking process. At this point, the calculation pro- 
ceeds via conservative molecular dynamics (MD) so as 
to allow hydrodynamic modes to develop. Here we use 
a time step 5t = 0. 005 ^/ma^Je. In the smallest system, 
crosslinking is carried out one bond at a time. A single 
gel particle is randomly selected and all other gel par- 
ticles within a distance of 1.2(To are identified. One of 
the particles in this list is randomly selected and bonded 
irreversibly to the central particle through the tethering 
potential Vnnirij) = ■^k{rij — tq)^ with k — 5e/aQ and 
rg = (7r/6^>)'^/'^CTo. Each gel particle is allowed to bond 
to no more than six others and bonding between any pair 
of particles occurs at most once. The configuration of the 
entire system is then updated for 100 time steps and the 
entire bonding process is repeated until ZpNgei crosslinks 
have been added. The parameter p is analogous to the 
occupation probability in a bond percolation process on 
the simple cubic lattice. In larger systems, the number 
of crosslinks added in the bonding steps is scaled by the 
system size in order to keep the crosslinking rate per gel 
particle constant. 

The parameters in the potentials and the total vol- 
ume fraction $ are the same as in our previous work 

. The differences are that in this earlier work all par- 
ticles were considered to be gel particles and that the 
crosslinking was done instantaneously, at i — 0, when 
the particles were on the vertices of a cubic lattice and 
thus all structural properties were those of percolation in 
three dimensions. The present model is similar in some 
ways to a model discussed by Gimel et al. |^ and Hasmy 
and JuUien | [lO| | who studied percolation in the context of 
diffusion-limited cluster-cluster aggregation using Monte 
Carlo methods. Their model differs from ours in that it 
is a lattice model, in the details of the crosslinking pro- 
cess, in the lack of solvent and in the nature of the cluster 
dynamics. In Monte Carlo simulations, one is forced to 
arbitrarily choose the mass-dependent diffusion constant 
D{m) whereas in our molecular dynamics calculations 
it is determined by the existing structure and the inter- 



particle forces. In the regime that is of interest here, i.e., 
high enough gel density that percolation is possible, these 
authors find the critical behavior of ordinary percolation. 

In a separate set of runs, we calculate the stress-stress 
autocorrelation function and, through the appropriate 
Green-Kubo formula, the shear viscosity. Equilibration 
and crosslinking are carried out as described above and 
the calculation of the viscosity is again done with conser- 
vative MD. 

The adjustable parameters in our calculations are the 
gel fraction c, the crosslink density p and the system size. 
Here we report results for c = 0.2, 0.3 and 1.0. Calcu- 
lations for other values of c are in progress and will be 
reported in a future publication |ll| . We parametrize the 
size of our system in terms of the dimensionless length 
L = iV^/'^ where N is the total number of gel and solvent 
particles. Because the crosslinking process is itself quite 
time consuming, we are able only to simulate systems up 
to size L — ?>2 (32768 particles) and this makes our es- 
timates of critical exponents rather imprecise. A second 
factor contributing to the uncertainty in critical expo- 
nents is that we need to determine the critical crosslink 
density Pc for each value of c whereas for lattice perco- 
lation this number is known to high accuracy. We next 
discuss the static (geometric) properties of our model. 



III. PERCOLATION 

The critical concentration pc at which percolation oc- 
curs in the thermodynamic limit L — > oo is accurately 
estimated from the intersection of curves /(p, L) as func- 
tion of p for different values of L. Here /(p, L) is the 
fraction of samples percolating in a system of size L at 
crosslink concentration p. For the two cases of interest 
here, c = 0.3 and c = 0.2, we find Pc = 0.3165 ± 0.0005 
and Pc = 0.3735 ± 0.001. Once pc has been determined, 
the correlation length exponent v can be estimated from 
the collapse of the data for the function / when plotted 
as function of {p — pc)L^^^ ■ We show this collapse of the 
data for c = 0.2 and 0.3 in figure (|l|). For c = 0.3, the best 
collapse of the data for 8 < L < 32 is obtained for v = 1.0 
which should be compared to the three-dimensional per- 
colation result V = 0.88. For c = 0.2, finite-size effects are 
more pronounced and the data for L = 8 have been ex- 
cluded. For this case, the best collapse of the data is ob- 
tained for u = 1.05. This method of estimating a critical 
exponent is not very accurate but the three-dimensional 
percolation value u — 0.88 provides a significantly worse 
collapse of the data. 

We next discuss the mean size of finite clusters be- 
cause this data provides an unbiased estimate of the ra- 
tio ^/v. In the thermodynamic limit, S{p) ~ \p — Pc\~"' 
with 7 w 1.8 for d = 3 percolation. For finite L, S{p, L) 
is peaked near pc with a peak height that grows as 
Therefore, rescaling the peak heights to the same value 
for different L provides an estimate of ^/v that is not 
affected by errors in either pc or v. Of course, the overall 
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FIG. 1: Fraction of samples percolating for c = 0.2 and c — 
0.3 as function of the scaled crosslink concentration x = {p ~ 
Pc)L}^'' for 8 < L < 32. The values of the exponent v used 
are 1.05 for c = 0.2 and 1.0 for c = 0.3. For L = 8 we 
have simulated 20000 independent crosslinkings at each p; for 
L = 32 the data are derived from 3000 samples for each p. 
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FIG. 2: Scaled form of the mean mass of finite clusters 
L-^/''S{p,L) for c = 0.3 and 8 < L < 32. Here -y/u = 1.815 
and I' = 1.0. 



collapse of the data to a universal curve depends on accu- 
rate determination of these two quantities but the peak 
height does not. In figures (||) and (||) we show the func- 
tion L~''/'^S{p, L) plotted as function of a; = {jp—pc)L^^^ 
for the previously determined values of Pc and v. The 
collapse to a universal curve is quite respectable for both 
c — 0.3 and 0.2 for 7/1/ = 1.815 and 1.80, respectively. 
As above, the data for L = 8 have been excluded for 
c = 0.2. We note that in the case of three-dimensional 
percolation the ratio ^jv k. 2.05. Use of this value of ^jv 
in figure {^) would result in a 40% difference between the 
peak heights for _L = 32 and L = 8. 

In the scaling theory of percolation 0], the ratio 7/1^ = 
(i(3 — ■'■)/(''" ~ 1), where d = 3 is the dimensionality and 
T is the exponent characterising the cluster size distri- 
bution at p = Pc- If we enforce this scaling relation, 
we obtain r « 2.25 for both c — 0.2 and c = 0.3. Us- 
ing a = [t — l)/di', we find <t{c = 0.3) — 0.415 and 
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FIG. 3: Same as figure in this case for c — 0.2 with 'yjv — 
1.80. 



ct(c = 0.2) = 0.417. Using D = II {av) for the fractal di- 
mension results in the prediction D{c = 0.3) = 2.41 and 
D{c = 0.2) = 2.29 for the fractal dimensions of the clus- 
ters. As well, the hyperscaling relation 2[3/v = 3 — 7/1^ 
yields iijv = 0.593 and 0.6 for c = 0.3 and 0.2 respec- 
tively. The accuracy of these scaling predictions is tested 
in figures (|) to (0). 

In figure we show the number of clusters n(m) of 
mass mz.tpK.pc for c ~ 0.2 and 0.3 for L = 32 and 
TO < 400. For the case m = 1, we have only counted 
the uncrosslinked gel particles. In neither case is the 
data well described by a simple power law, in contrast to 
percolation on a lattice or in the absence of solvent where 
the exponent r w 2.18 is already obtained for 2 < to < 
20. A fit to a power law over the range 20 < to < 400 
yields t = 2.13 for c = 0.2 and t = 2.16 for c = 0.3. The 
straight lines in figure (^) are the best fits to the form 
n = Am~'^-^^ over the range 20 < to < 400 and while the 
fit is not perfect, the data are not inconsistent with this 
behavior in the limit of large to. 

In figure we show the square of the radius of gyra- 
tion Rg (to) as function of to for a system of size L = 32 
together with curves to^/^^'^) with D{c = 0.2) = 2.29 
and D{c = 0.3) = 2.41 as determined above. The 
data again show considerable curvature but the fit to 
the assumed functional form is reasonable over the range 
20 < TO < 100. 

Finally, in figures (|6|) and we display the scaled 
form of P{L,p), the probability that a gel particle is part 
of the percolating cluster using the predicted exponent 
ratios P/i^ = 0.593 for c = 0.3 and = 0.6 for c = 0.2. 
These two figures present the least impressive collapse of 
data to a universal curve, especially at the larger values 
of P. One can improve the collapse by different choice of 
and V but at the expense of violating hyperscaling. 
We also note that the data for the two largest values of L 
are reasonably close to each other over the entire range 
of X. 

We have also carried out a limited number of simu- 
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FIG. 4: Number of clusters n(m) of mass m for p for c — 
0.2 and 0.3. The data for c = 0.2 has been lowered by a factor 
of 5 for separation of curves. The straight lines represent fits 
to am~^ with t determined by imposing hyperscaling (see 
text). 
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FIG. 5: Square of the radius of gyration Rg{m) as function 
of cluster mass m for p « Pc and c — 0.2 and 0.3. Straight 
Unes are fits to ]^{ni) — arr?^^ with the fractal dimensions 
determined by requiring that hyperscaling holds (see text). 



lations for c = 0.5 and c — 1.0 with the crosslinking 
process described above. In both cases, the critical ex- 
ponents and the cluster size distributions are entirely 
consistent with ordinary three-dimensional percolation. 
This suggests that either there is a critical gel fraction 
1^ below which the geometric properties of the clusters 
are described by continuously varying exponents or that 
the apparent variation of the exponents with c described 
above is a finite-size artefact. Only simulations of larger 
systems can resolve this issue. 



IV. SHEAR VISCOSITY 

We have calculated the shear viscosity for systems 
up to size L = 20 as function of the crosslink den- 



FIG. 6: Plot of the scaled form of the order parameter P{L,p) 
for pg = 0.3 and 8 < L < 32. The exponents are fijv = 0.593 
and V = 1.0. 
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FIG. 7: Same as figure (|) but for c = 0.2 and = 0.6 and 
V = 1.05. 



sity p for c = 0.3 and for L = 12 for c = 0.2. Sys- 
tems are equilibrated as a liquid, crosslinked as described 
above and then evolved by constant energy MD for 40000 
or 80000 time steps, depending on the crosslink den- 
sity. Here we have typically used 500 to 2000 differ- 
ent realizations of the crosslinks at each -p. We calcu- 
late, as in H, the stress-stress autocorrelation function 



'it) 



iEa</3(^"/3 (0^^0/3(0)) where 
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are elements of the stress tensor. Here the sum is over 
both gel and solvent particles and V is the derivative of 
the pair potential between particles i and j. The analy- 
sis of the stress-stress correlation function has been de- 
scribed in 1^ and is done in the same way here. As 
p Pc, Caa decays extremely slowly and is fitted, at 
long times, to a stretched exponential. The static shear 
viscosity is then obtained from the appropriate Green- 
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FIG. 8: The dimensionless shear viscosity 

aov{L,p)/{mkBTy^^ for c = 0.3 times L^"'^" plotted 
as function of x = {p — pc)L^^'^ with s/f = 0.425 and u = 1.0. 
The straight line represents the function a;""''" (see text). 
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FIG. 9: The dimensionless shear viscosity CTo»7(p)/(mfcsr)^''^ 
for L = 12 and c — 0.2 and 0.3 plotted as function of (pc — p). 
The straight line represents the function {pc — p)"" with s = 
0.425 for c = 0.3 and s = 0.3 for c = 0.2. 



Kubo formula [fl3| 



r] = lim 



1 



VkeT 



The results, for c = 0.3 are shown in finite-size scaled 
form in Fig. ^ where L'~ " / ri{L /£) is plotted as func- 
tion of the scaled concentration x . In contrast to our 
previous result for c — \ and instantaneous crosslinking 
where we found s « 0.7, we find that s k, 0.425 pro- 
vides an excellent collapse of the data with v = 1.0. We 
note that, outside the critical region, consistency of the 
finite-size scaling ansatz requires the scaled viscosity to 



vary as x 
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0.425 g_j-^jj jg (jiga,r that the data are 



consistent with this behavior. 

We have also calculated the shear viscosity for c = 0.2 
for L = 12. The raw data are displayed in figure (p|) as 



function of Pc—p together with the corresponding results 
for c = 0.3. Fitting to a power law outside the critical 
region produces an exponent s w 0.3 suggesting, as in 
the case of the static properties, a variation of critical 
exponents with c and an absence of universality. 



V. DISCUSSION 

In this article we have proposed and investigated a 
new model for gelation which incorporates a solvent on 
a microscopic level. For relatively small concentrations 
of gel, the geometric properties of the system close to 
the gel point seem to depend continuously on this gel 
fraction and are, at least for the system sizes investi- 
gated, markedly different from three-dimensional perco- 
lation. In particular, the fractal dimension of the clusters 
seems to be smaller than those of percolation clusters and 
this more spidery morphology may be responsible for the 
slower divergence of the shear viscosity as the gel point is 
approached. The change in the exponents controlling the 
geometric properties is rather small and further study of 
larger systems is certainly necessary to confirm this re- 
sult. However, the exponent s that characterizes the di- 
vergence of the shear viscosity at the gel point is reduced 
by almost a factor of 2 from its value in the absence of 
solvent and it is unlikely that this can be attributed to 
finite-size effects. In light of this result, it seems implau- 
sible that a single universality class describes the behav- 
ior of transport coefficients and, presumably, the moduli 
of the amorphous phase near the gel point. The con- 
siderable dispersion found in experimental values of the 
critical exponents IT^ is another indicator that this may 
be the case. 

In future work we intend to explore this new model in 
greater detail. It will be interesting to investigate if the 
exponent s and the static exponents are tunable by vary- 
ing the concentration of the solvent and the solubility of 
the solute. We also intend to study diffusion constants as 
function of cluster size and to investigate the existence of 
long time tails. Finally, one of the original motivations 
for this model is the existence of a body of experimental 
work that has yielded values in the range 1.1 — 1.3 for 
the viscosity exponent s. Clearly, we have moved further 
from this range of values compared to our previous re- 
sults. If the cluster size distribution and cluster geometry 
is the determining factor in the critical behavior of the 
transport coefficients then this indicates that models that 
produce more compact rather than more tenuous clusters 
than those arising from percolation may be appropriate. 

One of us (BJ) thanks the Physics Department at Si- 
mon Fraser University for its hospitality during a sab- 
batical visit. We thank Ralph Colby, Paul Goldbart and 
Sune Jespersen for helpful discussions. This research is 
supported by the NSERC of Canada. 



6 



[1] See for example D. Stauffer and A. Aharony, Introduction 
to Percolation Theory, 2nd Edition, (Taylor and Francis, 
London, 1994). 

[2] See D. StaufTer, A. Coniglio and M.Adam, Adv. Pol. Sci. 

44, 103 (1982) for a review. 
[3] D. Vernon, M. Plischke and B. Joos, Phys. Rev E 64, 

03105 (2001). 

[4] P.G. de Gennes, J. Phys (Paris), 40, L197 (1979). 

[5] K. Broderix, H. Lowe, P. Miiller and A. Zippelius, Eu- 

rophys. Lett., 48, 421 (1999); Phys. Rev. E 63, 011510 

(2001). 

[6] M. Adam, M. Delsanti, D. Durand, G. Hild and J.P. 
Munch, Pure Appl. Chem., 53, 1489 (1981); M. Adam, 
M. Delsanti and D. Durand, Macromolecules, 18, 2285 
(1985); D. Durand, M. Delsanti and J.M. Luck, Euro- 
phys. Lett., 3, 297 (1987). 

[7] CP. Lusignan, T.H. Mourey, J.C. Wilson and R.H. 
Colby, Phys. Rev. E 52, 6271 (1995); J.E. Martin and 
J. Wilcoxon, Phys. Rev. Lett., 61, 373 (1988); D. Adolf 
and J.E. Martin, Macromolecules, 23, 3700 (1990); J.E. 
Martin, J. Wilcoxon and J. Odinek, Phys. Rev. A 43, 
858 (1991); J.E. Martin, D. Adolf and J. Wilcoxon, Phys. 
Rev. Lett., 61, 2620 (1988). 

[8] J.P. Powles and D.M. Heyes, Mol. Phys. 98, 917 (2000). 
These authors have studied the properties of systems 
with a pair potential of the form Vifij) = e((T/rij)", 



for 12 < n < 288, including our case n = 36. 
[9] J.C. Gimel, D. Durand and T. Nicolai, Phys. Rev. B 51, 
11348 (1995). 

[10] A. Hasmy and R. Jullien, Phys. Rev. E 53, 1789 (1996). 

[11] M. Plischke, S. Jespersen and B. Joos, unpublished. 

[12] A critical volume fraction appears in the model of Gimel 
et al. [^], albeit in a rather trivial way: When the vol- 
ume fraction is increased above 0.31, the system perco- 
lates instantaneously (in the thermodynamic limit) and 
the features associated with cluster-cluster aggregation 
disappear. 

[13] M.P. Allen and D.J. Tildesley, Computer Simulation 
of Liquids (Oxford University Press, New York, 1987), 
Chap. 2; J.P. Hansen and I.R. MacDonald, Theory of 
Simple Liquids, 2nd ed. (Academic Press, New York, 
1986). 

[14] A fit of the raw data for i = 12 and p < 0.26 to the form 
ri{p) — a{pc—p)~'' produces the estimates pc = 0.310 and 
s = 0.37, in good agreement with the percolation theory 
estimate of pc and with the finite-size scaling analysis 
reported above. 

[15] For a review, see M. Adam and D. Lairez in The Physi- 
cal Properties of Polymeric Gels, J.P. Cohen Addad, ed. 
(John Wiley and Sons, Ltd., New York, 1996) p. 87. 



